#Cálculos estadísticos en hidrología
#Análisis de Frecuencia

#Lectura de datos 
#Importar datos desde Excel
#Libreria para cargar archivos excel
library(readxl)
#file.choose() 
filename<-"D:\\Paola\\2023B\\Hidrología Aplicada\\1. Análisis de Frecuencia\\Datos.xlsx"
max_preci<-read_excel(filename,skip = 1 )
summary(max_preci)
colnames(max_preci)<-c('año','d5','d10','d15','d20','d30','d60','d120','d360','d1440')

#Duracion a evaluar
preci_dmin<- max_preci$d1440

# Crear un Histograma Sencillo
hist(preci_dmin) 

#Histograma con librería ggplot
#Libreria para graficos
library(ggplot2)
ggplot(max_preci,aes(preci_dmin))+
  geom_histogram(bins=12, fill = "blue", color= 'white')+
  xlab('Rango de precipitacion [mm]')


# Distribución de Gumbel

#Paso 1 Obtencion de probabilidad de no ocurrencia ligada a un periodo de retorno
max_preci$p=1
max_preci$Tr=1
for(i in 1:nrow(max_preci)){
  max_preci[i,"p"]=1-(1/((nrow(max_preci)+1)/i))
  max_preci[i,"Tr"]= (nrow(max_preci)+1)/i
}

#Instalación de paquete para la utilizacion de distribución Gumbel
#install.packages("evd")
library(evd)

#Paso 2 Obtencion de los parámetros de la distribución
a<- mean(preci_dmin)-0.5463*(sd(preci_dmin)) #Parámetro de Gumbel
b<- sd(preci_dmin)/1.15 #Parámetro de Gumbel 

#Paso 3 Obtencion del cuantil que acumula la probabilidad -- Precipitacion Calculada
max_preci$pc<- qgumbel(max_preci$p,loc = a, scale = b)  


# Gráfico PObservada y Pcalculada

#Ordenar los datos de mayor a menor

max_preci$d1440or<- sort(preci_dmin, decreasing = TRUE)

#Grafico con ggplot
ggplot(max_preci,aes(x= max_preci$Tr, y= max_preci$d1440or))+ geom_point()+ 
  geom_line(aes(x = max_preci$Tr, y=max_preci$pc, col = "red"))+
  xlab("Periodo de retorno") + ylab("Precipitación")

#Ajuste y Prueba de bondad 
#Instalacion de paquetes 
#install.packages('fitdistrplus')
#install.packages("texmex")
library(fitdistrplus)
library(texmex)

#Ajuste a distribucion Gumbel
G <- fitdist(preci_dmin, "gumbel", 
                         start=list(mu=mean(preci_dmin), sigma=sd(preci_dmin)))
plot(G)
G$estimate

#Prueba de bondad Gumbel
PB_G<- gofstat(G)
PB_G$ks
PB_G$chisqpvalue


#Resumen de bondad de los ajustes

gofstat(G,fitnames = "Gumbel")

#Distribucion de Gumbel paquete texmex

qgumbel(max_preci$p,mu = mean(preci_dmin),sigma = sd(preci_dmin))
